Hydrodynamic synchronisation of non-linear oscillators at low Reynolds number. 
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We introduce a generic model of weakly non-linear self-sustained oscillator as a simplified tool 
to study synchronisation in a fluid at low Reynolds number. By averaging over the fast degrees of 
freedom, we examine the effect of hydrodynamic interactions on the slow dynamics of two oscillators 
and show that they can lead to synchronisation. Furthermore, we find that synchronisation is 
strongly enhanced when the oscillators are non-isochronous, which on the limit cycle means the 
oscillations have an amplitude-dependent frequency. Non-isochronity is determined by a nonlinear 
coupling a being non-zero. We find that its (a) sign determines if they synchronise in- or anti- 
phase. We then study an infinite array of oscillators in the long wavelength limit, in presence 
of noise. For a > 0, hydrodynamic interactions can lead to a homogeneous synchronised state. 
Numerical simulations for a finite number of oscillators confirm this and, when a < 0, show the 
propagation of waves, reminiscent of metachronal coordination. 



Collections of cilia and flagella are examples of systems 
that display synchronisation [![. They are microscopic 
active filaments attached to the membrane of pro- and 
eukaryote cells Q whose synchronisation is thought to aid 
the efficiency of transport at the cellular scale. Typically 
arrays of cilia generate fluid flows along tissues but can 
also be used, like flagella, for the self-propulsion of swim- 
ming cells. Due to their tiny size, the Reynolds number 
associated with these flows is negligible. The coordinated 
beating of cilia is also thought to have important devel- 
opmental implications, such as the left-right symmetry 
breaking in the arrangement of the internal organs in the 
early embryo Q . A precise understanding of the role hy- 
drodynamics plays in their synchronised motion, is still 
missing. 

Both cilium and flagellum are made of complex sub- 
units, microtubules driven by molecular motors, and 
their modelling can be done at many levels. As syn- 
chronisation takes place on length-scales larger than the 
individual filaments, to a first approximation the fine 
details of their internal structure can be ignored. This 
coarse-grained approach has led to model studies of self- 
sustained oscillators rotating beads 043; beating 
filaments Q, as well as rigid rotating helices, [1, fl"o| . 
More recent work has focused on the conditions for hy- 
drodynamic synchronisations for two oscillators and 
the phase dynamics of oscillators with long range inter- 
Relatcd experiments investigating the dy- 
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namics of micro-systems have been performed in vivo on 
algae, [l3|,[l4| and on simple model systems [lj|, and even 
a macroscopic scale model of rotating paddles [r| . All 
these studies suggest that simple forms of active forces, 
e.g. as prescribed functions of time, are not enough 
to guarantee synchronisation. Rather, a complex, non- 
linear relation between forces and velocities is necessary. 
Important questions therefore are what aspects of hydro- 
dynamic interactions aid synchronisation and what fea- 
tures of oscillators make them good hydrodynamic syn- 
chronizers. 

The dynamics of a system close to an oscillatory in- 



stability can be conveniently described by weakly non- 
linear oscillators whose averaged equations are univer- 
sal [l| . This implies that the long time behaviour of many 
systems with simple spontaneous oscillations can be cap- 
tured by a generic model with a few parameters. Using 
this insight, in this paper we introduce a minimal model 
of an oscillator at low Reynolds number. To simplify our 
presentation, we study our model in one-dimension. At a 
coarse grained level, this degree of freedom can be inter- 
preted as the centre of a filament beating in a plane 17 1 . 

The slow dynamics of the oscillator is naturally char- 
acterised using of two variables: the amplitude and the 
phase. Under arbitrary initial conditions, the trajecto- 
ries of an isolated oscillator on long timescales converge 
to a closed curve, the limit cycle [lj|- While the ampli- 
tude is tightly constrained to the limit cycle curve, the 
phase can vary more freely. Hence many model studies 
of synchronisation have focused only on the phase dy- 
namics EH G3 • Our goal in this paper is to anal- 
yse the role played by both the amplitude and phase dy- 
namics on phase synchronisation mediated by hydrody- 
namics. We first study a pair of well separated deter- 
ministic oscillators and find that hydrodynamic interac- 
tions strongly enhance phase locking, if the oscillations 
are non-isochronous, which on the limit cycle means that 
the frequency of oscillations depends on the amplitude. 
We then consider an array of many oscillators, still well 
separated, in the presence of fluctuations. On long wave- 
lengths their slow dynamics can be naturally represented 
in terms of a broken symmetry (phase) variable, which 
is a non-equilibrium analogue of a Goldstone mode [l9| . 
Denoting by a ^ the parameter responsible for the non- 
isochronity of the oscillations, we find that when a > 0, 
hydrodynamic interactions can lead to in-phase synchro- 
nisation of the array. These results are confirmed by 
numerical simulations, which show also that conversely, 
for a < 0, the synchronisation is more subtle and leads 
to the propagation of waves. 

a. The model oscillator A universal model for stable 
spontaneous oscillations is provided by the normal form 
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of a dynamical system close to a supercritical Hopf bifur- 
cation [IH . To be concrete, we represent the oscillator in 
a low Reynolds number fluid as a sphere of radius a sub- 
ject to a time-varying force /. The equation of motion 
for the sphere, with x its deviation from its equilibrium 
position, is 
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where 7 = 67r?7a is the Stokes drag. The dynamics is 
encoded in the evolution equation for the force / : 
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Here, all the parameters, except a are positive quantities, 
The 1st and 3rd term of eq ([2]) give rise to respectively, 
a linear and a non-linear passive oscillator, while the 2nd 
term is responsible for active, self-sustained oscillations. 
We emphasize that all the terms in eq would emerge 
naturally from coarse-graining any friction-dominated 
microscopic model oscillator [J, [13, [13 • Eqs (TTJ), ^ 
can be conveniently non-dimcnsionaliscd as x = /; and 



FIG. 1: (color online) One dimensional lattice of oscillators. 
The inset illustrates the dynamic variables of a pair. 



To proceed, we take the time derivative of both sides 
of eq ([3]) and use, on the rhs, the evolution equation for 
the forces and the expression of forces as functions of 
velocities Xi obtained by inverting eq ([3]) as an expan- 
sion in a/r. Thus, to leading order, wc obtain equations 
for oscillators with reactive couplings [l[ (given by ® ) 

as i\ + lo 2 x\ 



6^/(1 - x 2 ) + e a x 3 , choosing units where by u; 2 = — . Note terms like -^^(^,7^ 



/ = - 

r = ?. They correspond to a weakly non- linear Van 
der Pol-Duffing oscillator [18j . The parameters := j 
and e a := |J are small quantities. We restrict ourselves 
here mainly to cases where e a / 1 = 0{1). 

b. Two oscillators coupled hydrodynamically The 
oscillators are arranged along the x-axis. The forces fi 
acting on the spheres, for i = 1,2, arc directed along 
the same axis and cause sphere 1 to oscillate around the 
origin and sphere 2 around position d. We denote by x% 
the deviations from these equilibrium positions, see fig[T] 
Their equations of motion are 



±J 7 1 (x 1 ,jx 1 ) + ^x 2 and x 2 + 0JqX 2 = 
— J r 2(x2, 7x2) + where D := —H(d)K loq represents 
the natural frequency of the linear oscillators, defined 
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xx = i (/1 + H{r)f 2 ) ■ 
i2 = ±(h + H(r)f 1 ), 



(3) 



C(edep) and ^f/i, of order O(e^) have been neglected. 

We now derive the equations governing the slow dy- 
namics of the oscillators [l|. This is done naturally 
using a complex amplitude A k and its complex con- 
jugate A k related to position and velocities by Xk = 
\{A k e iult + A* k e~ luJt ) and x k = l -f{A k e luJt - A^e - *"*) for 
k = 1,2. This requires of course that A* k = -A k e 2tuJt . 
Here ui is the (unknown) frequency of the non-linear os- 
cillators, determining the period, T = of the (fast) os- 
cillations. The (slow) complex amplitudes, on the other 
hand, hardly change on this timescale. Writing eq ^ 
and the dynamic equations for the forces in terms of A k 
and A* k and averaging over the period T one obtains 

Ai = -iAAi + AAi - (/3 + ix)A 1 \A 1 \ 2 + iSA 2 



where H (r) is a scalar, representing the hydrodynamic in- 
teractions, and r := d + x 2 —x\ is the separation between 
the sphere centres. We shall consider the limit of large 
separation r compared to the sphere radius a. Then, for 
an unbounded three-dimensional fluid, interactions are 
described by the Oseen tensor [2(| as H(r) = I 2 . For a 
rigid surface with a non-slip boundary condition, placed 
at distance h from the oscillators, one obtains effective in- 
teractions scaling as H(r) ~ [2l|. For an assembly of 
oscillators arranged on a regular lattice, d can be thought 
of as the lattice spacing, see fig [1] We assume that it 
is large compared to the amplitude of the oscillations, 
d > X2 — xi, and that the ratio ed '■= a/d, characterising 
the hydrodynamic coupling, satisfies td -C e M , t a - The 
time evolution of forces is given by fi = ^(fi,Xi), with 
ty(fi,Xi) defined in eq ([2]), and is entirely local j4j]. The 
long-range hydrodynamic coupling links the coordinates 
Xi via eq ([3]). In the following we denote the nonlinear 
parts of V(fi,Xi) by T^x^fi) := ^/ f (l - ax 2 ) + ax\. 
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Writing the complex amplitudes A k in polar form, 
A k = R k e i ^' k , eqs (j4| become a coupled system for the 
amplitudes R k and the phases <p k . Finally, this system 
can be reduced to a single equation for the phase dif- 
ference [l|. This can be achieved perturbatively, when 
the parameter <5, parametrising the hydrodynamic inter- 
actions, is small compared to the other terms. If in- 
teractions are neglected, R k have fixed points given by 

Rk -- 



■| . The dynamics of small deviations from these 



s k 



fixed points can be studied by writing R k - 

for st « 1. One finds that the deviations s k relax quickly 
to zero. Setting s k = we obtain s k as functions of 
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the phase difference tp := <p2 — (pi- The resulting expres- 
sions are then substituted in the equations for the phases. 
From them one obtains an Adler equation [l[ for tp, 

tp = v — 2— sin^/>. (5) 

Hence, cq ((5]) illustrates that phase locking is deter- 
mined by the hydrodynamic coupling, via 5, provided 
the oscillator is nonisochronous, i.e. a^O. Note that 
scales as ~ - — ed and v is related to the difference of the 
natural frequencies of the oscillators. We choose them to 
be identical, so we can set v = 0. While for v ^ varying 
the ratio of v and ^ controls the saddle-node bifurcation 
of cycles [lU, for v = eq (J5J) has a stable fixed point 
given by one of the zeros of sin(^) for tp G [0, 2tt]. The 
position of the stable point is determined by the sign of 
— which in turn is determined solely by the sign of 
the non-isochronism parameter a: when a < 0, then the 
equation has a stable fixed point at tp = it, i.e. the os- 
cillators lock in anti-phase; vice- versa, if a > then the 
equation has a stable fixed point at tp = and the os- 
cillators lock in-phase. A numerical solution, using the 
Euler method, of eq ([3]) confirms this. 

It is also interesting to note that the two fiagella of 
the microscopic algae C. Reinhardtii are found to alter- 
nate between periods of synchronised (with small phase 
difference) and non-synchronized beating [l]| [l4[ • This 
is well described by a stochastic Adler equation, of the 
same form as eq ([5]) but with an additional fluctuating 
term [3]. The estimates of the parameters presented 
in [l4| , for the flagellar synchronisation, indicate positive 
values for a and v of our model. 

When a = 0, we need to include higher order correc- 
tions in deriving eq ([5]) . Upon doing this we find to lead- 
ing order tp as —Sede^l+j^-costp] s'mip. When e d < 

the synchronisation is in-phase. Otherwise, both in- and 
anti-phase states are possible and synchronization de- 
pends on details such as initial conditions (confirmed nu- 
merically). These higher order terms also indicate that 
the transition from in-phase to anti-phase in general oc- 
curs at some a c ^ 0. Unsurprisingly when a = 0, syn- 
chronisation occurs more slowly (a higher order effect). 

c. Many oscillators coupled hydrodynamically As we 
have discussed above, the amplitudes of the oscillators 
are tightly constrained to the limit cycle and the long 
time behaviour can be reduced to an effective (amplitude 
dependent) dynamics of the phases. For a large num- 
ber N of oscillators, in the dilute regime, this is done 
by introducing the one-particle probability c(ip, y, t) = 
(fEL S(.V-<t>k{t))5(y-y k (t))) of having an oscillator 
with slow phase (p, at site y at time t, where the brack- 
ets () indicate the average over noise. The probability 
satisfies a Smoluchowski equation 

d tC = Dd 2 w c-d v {[u l+ n]c). (6) 

D is the diffusion coefficient resulting from both thermal 
and active fluctuations, wi the deterministic contribution 



of an isolated oscillator with ui\ = — A — and fi the 
deterministic effect of the hydrodynamic interactions, 

tt{y,V,t) ■= J dy 2 dip2c(ip 2 ,y2,t)<p mt (y2-y,<p,V2)- (7) 

<p mt = 2^- sm.{f2 — <p) + S' cos((^2 — <p) is obtained from 
the dynamics of two oscillators, (see eq (J5J). It describes 
the effect of the interactions on the phase of one oscillator 
due to the presence of the another. Here, 8' := 8(\y2 — y\)- 
The 1-particle probability can be expressed as an ex- 
pansion in its moments: 

<<P, y>t) = ^ [P(f> *) + (e^<%, t) + c.c.) + . . .] (8) 
To study synchronization we only need the first two : 

p(y,t) := dtpc{tp,y,t) ; (density) 
Jo 

$>(y,t) := / dipe lv c(ip,y,t) ; (1st harmonic) . (9) 
Jo 

The emergence (or not) of a globally synchronized state 
is obtained from the homogeneous probability c°(ip,t), 
with associated moments p°(t), representing spa- 

tially homogeneous dynamical states. The correspond- 
ing expression for O is obtained by evaluating the 
space integral in eq ([7]) with c = c . For hydrody- 
namic interactions scaling as H{r) ~ ^ the leading term 
from the integral depends both on the lattice spacing 
d, and the total length L of the array. Hence, £l°(t) = 
^-ln(i/d)[-i| + l]e _l¥,1 $°(i) + c.c. For interactions 

scaling as H(r) ~ the leading term in the integral 
depends only on the lattice spacing d. Consequently, the 
term a In (L/d) is replaced by one ~ 9 ^ r . Dynamic equa- 
tions for the homogeneous moments are derived by tak- 
ing the time derivative of both sides of eq © , inserting 
cq ((6]) and using eq dHJ to close the system. Since p is a 
conserved variable, dtp = 0, while $° satisfies 

a t $° = r$°. (10) 

It is worth noting that in the absence of noise 
c°(<p,t) = ^£f =1 <5(<^ - fk{t)) and $°(t) reduces to 
the order parameter introduced by Kuramoto, 3>°(t) = 
■h £feLx e l<Pk ^ t \ representing the (mean field) average 
over a population of oscillators [l|, [TU, [22[ . 

It is useful to express $°(f) = P°(t)e lQ0< -^ in polar 
form (reflecting the U(l) symmetry). We obtain equa- 
tions for its amplitude and phase as dtP° = Re[T]P° 
and d t Q° = Im[T}. Re(F) = -(£> - f ^ \ n (L/d)p°) 
is the real part of T. Here, the first term is due to 
noise, whereas the second term encodes the effect of 
two body interactions. The imaginary part is Im(T) = 
\^ + ^HL/d)p Q ]. 

As in the Kuramoto model 0, [22| , order (synchronisa- 
tion) is determined by a non-zero, constant value of P . 
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FIG. 2: (color online) Space-time plots of the positions for 
N — 100 deterministic oscillators, (D — 0). (a), (b) describe 
respectively the case for a > and a < 0, after long time. 
The initial conditions of the oscillators are the same for both 
values of a: identical amplitudes, close to the maximum value, 
and random, Gaussian distributed, phases. The parameters of 
the model are 7 = 10" 3 Pa s urn; * = l J ~; u = 0.05^; 
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and a/d « 0.005. 



Here, the dynamic equation for P° shows that the onset 
of order is controlled by the sign of i?e[r]. If Re[T] < 0, 
order is suppressed. On the contrary, when i?e[r] > 0, or- 
der is enhanced. A stabilising term of the type ~ $ |$°| 2 
in eq (fTU|) is needed for P° to stop unbounded growth and 
attain a finite value at long times. Such a term could be 
generated for instance by taking into account three-body 
interactions. Finally, the condition Re[T] = defines 
a transition line in the space of parameters [l9| . Cru- 
cially, from these considerations, homogeneous synchro- 
nization is possible only when a > 0: (i) in presence of 
noise (D ^ 0) and by keeping all the parameters fixed, 
synchronisation occurs only above a particular value of 
density; (ii) neglecting noise (D = 0), instead, synchroni- 
sation occurs for any (finite) value of the density. On the 
contrary, when a < both terms in Re(T) are negative 
and homogeneous order is prohibited. This behaviour 
suggests a spin analogy, where a > (ferromagnet) pro- 
motes alignment of neighbouring oscillator phase (spins) 



while a < (antifcrromagnct) promotes anti-alignment. 

We compared these results with numerical simulation 
for a large but finite number of deterministic oscillators 
{D = 0). In fig [2] we show typical space-time plots for the 
positions of N = 100 oscillators and compare the effects 
of different signs of a. For a > 0, sec fig[2ja), the sys- 
tem displays spatially homogeneous order, i.e. in-phase 
synchronised state. Interestingly, when a < 0, although 
homogeneous order is lacking, fig EJb) still shows a co- 
herent motion of the oscillators, with propagating waves. 
As suggested by the antifcrromagnctic analogy, the oscil- 
lators self-organise into a dynamical state which is close 
to the anti-phase synchronised state, but deviates from 
it at long wavelengths. 

In conclusion, we have presented a simple, one- 
dimensional model (that can be generalised to higher di- 
mensions [23|) to investigate analytically the role of hy- 
drodynamic interactions on the synchronisation dynam- 
ics of oscillators at low Reynolds number. We studied the 
case of two oscillators and found that synchronisation, ei- 
ther in- or anti-phase, was determined to leading order by 
both hydrodynamic interactions and non-isochronism of 
the oscillations (a ^ 0). We then derived a coarse grained 
description for an infinite array of oscillators and found 
that spatially homogeneous order, corresponding to the 
in-phase synchronisation of the array, can occur only for 
a > 0. Systems of cilia are known to display metachronal 
waves [24| . Our analysis suggests that these could be ob- 
tained in two different ways: either as slow hydrodynamic 
(phase) modes, like spin waves, when a > 0; or alterna- 
tively, for a < 0, as a spatially inhomogeneous, approxi- 
mately anti-phase synchronised state, as indicated by the 
numerics. A more extensive investigation of these issues 
is left for the future. 
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